home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / perlin.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  8KB  |  226 lines
  1. /*
  2.  *    Noise and Dnoise routines
  3.  *
  4.  *    Many thanks to Jon Buller of Colorado State University for these
  5.  *    routines.
  6.  */
  7. #include "vect.h"
  8.  
  9. #define NUMPTS    512
  10. #define P1    173
  11. #define P2    263
  12. #define P3    337
  13. #define phi    0.6180339
  14.  
  15. static double pts[NUMPTS];
  16. static firsted;
  17.  
  18. static initnoise()
  19. {
  20.    int i;
  21.    
  22.    if(!firsted) {
  23.            for (i = 0; i < NUMPTS; ++i)
  24.       pts[i] = (random()%1000) / 1000.0;
  25.     firsted = 1;
  26.     }
  27. }
  28.  
  29. double pnoise(p)
  30. vect *p;
  31. {
  32.    int xi, yi, zi;
  33.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  34.    double xf, yf, zf;
  35.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  36.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  37.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  38.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  39.  
  40.    initnoise();
  41.    xi = floor(p->x);
  42.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  43.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  44.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  45.  
  46.    yi = floor(p->y);
  47.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  48.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  49.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  50.  
  51.    zi = floor(p->z);
  52.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  53.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  54.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  55.  
  56.    p000 = pts[xa + ya + za & NUMPTS - 1];
  57.    p100 = pts[xb + ya + za & NUMPTS - 1];
  58.    p200 = pts[xc + ya + za & NUMPTS - 1];
  59.    p010 = pts[xa + yb + za & NUMPTS - 1];
  60.    p110 = pts[xb + yb + za & NUMPTS - 1];
  61.    p210 = pts[xc + yb + za & NUMPTS - 1];
  62.    p020 = pts[xa + yc + za & NUMPTS - 1];
  63.    p120 = pts[xb + yc + za & NUMPTS - 1];
  64.    p220 = pts[xc + yc + za & NUMPTS - 1];
  65.    p001 = pts[xa + ya + zb & NUMPTS - 1];
  66.    p101 = pts[xb + ya + zb & NUMPTS - 1];
  67.    p201 = pts[xc + ya + zb & NUMPTS - 1];
  68.    p011 = pts[xa + yb + zb & NUMPTS - 1];
  69.    p111 = pts[xb + yb + zb & NUMPTS - 1];
  70.    p211 = pts[xc + yb + zb & NUMPTS - 1];
  71.    p021 = pts[xa + yc + zb & NUMPTS - 1];
  72.    p121 = pts[xb + yc + zb & NUMPTS - 1];
  73.    p221 = pts[xc + yc + zb & NUMPTS - 1];
  74.    p002 = pts[xa + ya + zc & NUMPTS - 1];
  75.    p102 = pts[xb + ya + zc & NUMPTS - 1];
  76.    p202 = pts[xc + ya + zc & NUMPTS - 1];
  77.    p012 = pts[xa + yb + zc & NUMPTS - 1];
  78.    p112 = pts[xb + yb + zc & NUMPTS - 1];
  79.    p212 = pts[xc + yb + zc & NUMPTS - 1];
  80.    p022 = pts[xa + yc + zc & NUMPTS - 1];
  81.    p122 = pts[xb + yc + zc & NUMPTS - 1];
  82.    p222 = pts[xc + yc + zc & NUMPTS - 1];
  83.  
  84.    xf = p->x - xi;
  85.    x1 = xf * xf;
  86.    x2 = 0.5 * x1;
  87.    x1 = 0.5 + xf - x1;
  88.    x0 = 0.5 - xf + x2;
  89.  
  90.    yf = p->y - yi;
  91.    y1 = yf * yf;
  92.    y2 = 0.5 * y1;
  93.    y1 = 0.5 + yf - y1;
  94.    y0 = 0.5 - yf + y2;
  95.  
  96.    zf = p->z - zi;
  97.    z1 = zf * zf;
  98.    z2 = 0.5 * z1;
  99.    z1 = 0.5 + zf - z1;
  100.    z0 = 0.5 - zf + z2;
  101.  
  102.    return   z0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  103.                   y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  104.                   y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  105.             z1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  106.                   y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  107.                   y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  108.             z2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  109.                   y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  110.                   y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  111. }
  112.  
  113. Dpnoise(p,dn)
  114. vect *p, *dn;
  115. {
  116.    int xi, yi, zi;
  117.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  118.    double xf, yf, zf;
  119.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  120.    double xd2, xd1, xd0, yd2, yd1, yd0, zd2, zd1, zd0;
  121.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  122.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  123.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  124.  
  125.    initnoise();
  126.    xi = floor(p->x);
  127.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  128.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  129.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  130.  
  131.    yi = floor(p->y);
  132.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  133.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  134.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  135.  
  136.    zi = floor(p->z);
  137.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  138.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  139.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  140.  
  141.    p000 = pts[xa + ya + za & NUMPTS - 1];
  142.    p100 = pts[xb + ya + za & NUMPTS - 1];
  143.    p200 = pts[xc + ya + za & NUMPTS - 1];
  144.    p010 = pts[xa + yb + za & NUMPTS - 1];
  145.    p110 = pts[xb + yb + za & NUMPTS - 1];
  146.    p210 = pts[xc + yb + za & NUMPTS - 1];
  147.    p020 = pts[xa + yc + za & NUMPTS - 1];
  148.    p120 = pts[xb + yc + za & NUMPTS - 1];
  149.    p220 = pts[xc + yc + za & NUMPTS - 1];
  150.    p001 = pts[xa + ya + zb & NUMPTS - 1];
  151.    p101 = pts[xb + ya + zb & NUMPTS - 1];
  152.    p201 = pts[xc + ya + zb & NUMPTS - 1];
  153.    p011 = pts[xa + yb + zb & NUMPTS - 1];
  154.    p111 = pts[xb + yb + zb & NUMPTS - 1];
  155.    p211 = pts[xc + yb + zb & NUMPTS - 1];
  156.    p021 = pts[xa + yc + zb & NUMPTS - 1];
  157.    p121 = pts[xb + yc + zb & NUMPTS - 1];
  158.    p221 = pts[xc + yc + zb & NUMPTS - 1];
  159.    p002 = pts[xa + ya + zc & NUMPTS - 1];
  160.    p102 = pts[xb + ya + zc & NUMPTS - 1];
  161.    p202 = pts[xc + ya + zc & NUMPTS - 1];
  162.    p012 = pts[xa + yb + zc & NUMPTS - 1];
  163.    p112 = pts[xb + yb + zc & NUMPTS - 1];
  164.    p212 = pts[xc + yb + zc & NUMPTS - 1];
  165.    p022 = pts[xa + yc + zc & NUMPTS - 1];
  166.    p122 = pts[xb + yc + zc & NUMPTS - 1];
  167.    p222 = pts[xc + yc + zc & NUMPTS - 1];
  168.  
  169.    xf = p->x - xi;
  170.    x1 = xf * xf;
  171.    x2 = 0.5 * x1;
  172.    x1 = 0.5 + xf - x1;
  173.    x0 = 0.5 - xf + x2;
  174.    xd2 = xf;
  175.    xd1 = 1.0 - xf - xf;
  176.    xd0 = xf - 1.0;
  177.  
  178.    yf = p->y - yi;
  179.    y1 = yf * yf;
  180.    y2 = 0.5 * y1;
  181.    y1 = 0.5 + yf - y1;
  182.    y0 = 0.5 - yf + y2;
  183.    yd2 = yf;
  184.    yd1 = 1.0 - yf - yf;
  185.    yd0 = yf - 1.0;
  186.  
  187.    zf = p->z - zi;
  188.    z1 = zf * zf;
  189.    z2 = 0.5 * z1;
  190.    z1 = 0.5 + zf - z1;
  191.    z0 = 0.5 - zf + z2;
  192.    zd2 = zf;
  193.    zd1 = 1.0 - zf - zf;
  194.    zd0 = zf - 1.0;
  195.  
  196.    dn->x      = z0 * (y0 * (xd0 * p000 + xd1 * p100 + xd2 * p200) +
  197.                       y1 * (xd0 * p010 + xd1 * p110 + xd2 * p210) +
  198.                       y2 * (xd0 * p020 + xd1 * p120 + xd2 * p220)) +
  199.                 z1 * (y0 * (xd0 * p001 + xd1 * p101 + xd2 * p201) +
  200.                       y1 * (xd0 * p011 + xd1 * p111 + xd2 * p211) +
  201.                       y2 * (xd0 * p021 + xd1 * p121 + xd2 * p221)) +
  202.                 z2 * (y0 * (xd0 * p002 + xd1 * p102 + xd2 * p202) +
  203.                       y1 * (xd0 * p012 + xd1 * p112 + xd2 * p212) +
  204.                       y2 * (xd0 * p022 + xd1 * p122 + xd2 * p222));
  205.                                   
  206.    dn->y      = z0 * (yd0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  207.                       yd1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  208.                       yd2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  209.                 z1 * (yd0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  210.                       yd1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  211.                       yd2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  212.                 z2 * (yd0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  213.                       yd1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  214.                       yd2 * (x0 * p022 + x1 * p122 + x2 * p222));
  215.                                   
  216.    dn->z      = zd0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  217.                        y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  218.                        y2 * (x0 * p020 + x1 * p120 + x2 * p220)) +
  219.                 zd1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  220.                        y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  221.                        y2 * (x0 * p021 + x1 * p121 + x2 * p221)) +
  222.                 zd2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  223.                        y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  224.                        y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  225. }
  226.